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Recently exploratory studies were performed on the possibility of constraining the neutron star 
equation of state (EOS) using signals from coalescing binary neutron stars, or neutron star-black hole 
systems, as they will be seen in upcoming advanced gravitational wave detectors such as Advanced 
LIGO and Advanced Virgo. In particular, it was estimated to what extent the combined information 
from multiple detections would enable one to distinguish between different equations of state through 
hypothesis ranking or parameter estimation. Under the assumption of zero neutron star spins both 
in signals and in template waveforms and considering tidal effects to 1 post-Newtonian (IPN) order, 
it was found that 0(20) sources would suffice to distinguish between a stiff, moderate, and soft 
equation of state. Here we revisit these results, this time including neutron star tidal effects to 
the highest order currently known, termination of gravitational waveforms at the contact frequency, 
neutron star spins, and the resulting quadrupole-monopole interaction. We also take the masses 
of neutron stars in simulated sources to be distributed according to a relatively strongly peaked 
Gaussian, as hinted at by observations, but without assuming that the data analyst will necessarily 
have accurate knowledge of this distribution for use as a mass prior. We find that especially the 
effect of the latter is dramatic, necessitating many more detections to distinguish between different 
EOS and causing systematic biases in parameter estimation, on top of biases due to imperfect 
understanding of the signal model pointed out in earlier work. This would get mitigated if reliable 
prior information about the mass distribution could be folded into the analyses. 

PACS numbers: 26.60.Kp, 95.85.Sz 


I. INTRODUCTION 

Second-generation ground-based interferometric grav¬ 
itational wave (GW) detectors are currently under con¬ 
struction: Advanced LIGO [1] in the US, Advanced Virgo 
[5] in Italy, and KAGRA [3] in Japan. GEO-HF in Ger¬ 
many nig is already taking data. Later in the decade, 
LIGO-India [5] may join this global network of obser¬ 
vatories. Among the most promising sources for a first 
direct detection of gravitational waves are compact bi¬ 
naries composed of neutron stars or black holes, with 
detection rates in the range 1 — 100 yr“^ depending on 
the astrophysical event rate, the instruments’ duty cycle, 
and the sensitivities of the detectors dll]; see also [3] for 
detection rates under the assumption that short, hard 
gamma ray bursts are caused by coalescing binaries. 

Coalescing binaries consisting of two neutron stars 
(BNS), a neutron star and a black hole (NSBH), or two 
black holes (BBH) have rich scientific potential. They 
can be used to test general relativity in the genuinely 
strong-held regime [inHH] and are self-calibrating “stan- 
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dard sirens” for cosmology [TSUIK] . Moreover, BNS and 
NSBH coalescences can be used as probes of the elusive 
neutron star equation of state (EOS), about which little 
is currently known m- 

The possibility of constraining or measuring the neu¬ 
tron star EOS with gravitational wave observations of 
BNS coalescences has recently been the subject of ex¬ 
tensive investigation. The way in which the EOS en¬ 
ters the GW signal from a coalescing binary is mainly 
through tidal deformation. During the last stages of 
inspiral, the tidal held £ij of one component of a bi¬ 
nary will induce a quadrupole moment Qij in the other, 
where to leading order in the adiabatic approxima¬ 
tion Qij = —X(EOS;m) £ij. The tidal deformability 
A(EOS;to) depends on the neutron star mass m in a 
way that is governed by the EOS. This deformation 
of the neutron stars has an effect on the orbital mo¬ 
tion, and hence on the waveform of the emitted grav¬ 
itational wave signal; in particular, it enters the phase 
$(t). The deformability is related to the radius i?(m) 
through A(m) = {2/3)k2{m)R^{m), where ^2 is the sec¬ 
ond Love number. Although tidal effects enter the phase 
at high apparent post-Newtonian order (first appear¬ 
ing alongside the 5 post-Newtonian (5PN) phase con¬ 
tribution), these corrections come with a large prefactor: 
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\{m)/M^ cx (R/M)^ ~ 10^ — 10^ (with M the total mass 
of the binary), so that they may be observable with ad¬ 
vanced detectors. We stress that there are other ways 
in which the EOS enters the gravitational waveform; we 
will come to these momentarily. 

Read et al. estimated that with a single close-by source 
(at a distance of ^ 100 Mpc), the neutron star radius 
could be constrained to 10 % m Hinderer et al. per¬ 
formed a Fisher matrix analysis with PN waveforms trun¬ 
cated at 450 Hz to see how well the neutron star tidal de- 
formability might be measurable from the low-frequency 
inspiral dynamics alone |18j , concluding that it would be 
difficult to extract much information from this regime, at 
least with second-generation detectors. Damour, Nagar, 
and Villain performed a Fisher matrix calculation using 
effective one-body waveforms up to the point where the 
neutron stars are touching each other [19], which sug¬ 
gested that it might be possible after all to gain informa¬ 
tion about the EOS with advanced detectors. Lackey et 
al. performed a similar analysis for NSBH, which also 
indicated encouraging prospects |2Q|. The abovemen- 
tioned work considered single detections; a first study 
for multiple detected sources was performed by Markakis 
et al. | 21 | . who concluded that a similar accuracy to m 
could be attained with 3 sources that have low signal- 
to-noise ratio (SNR). On the other hand, Fisher matrix 
estimates can be unreliable at low SNR , prompt¬ 

ing more in-depth assessments. 

The first fully Bayesian investigation of the problem, 
in a realistic data analysis setting, was performed by Del 
Pozzo et al. [53]. BNS signals were “injected” into sim¬ 
ulated detector noise, assuming the projected final de¬ 
sign sensitivity of the Advanced LIGO-Virgo detector 
network. Sources were distributed in an astrophysically 
realistic way, leading to the distribution of SNRs that we 
expect to see towards 2018. Two different Bayesian anal¬ 
ysis methods were employed: hypothesis ranking within 
a list of different theoretical EOS, and parameter estima¬ 
tion. The former trivially allows one to combine infor¬ 
mation from multiple sources so as to arrive at a stronger 
result. To do the same with the latter method, parame¬ 
ters need to be identified that do not vary from source to 
source; in these were simply taken to be coefficients 
in a Taylor expansion of X{m) in powers of (m —mo)/M 0 , 
where uiq is some reference mass. A similar analysis to 
[25| in terms of parameter estimation was recently per¬ 
formed by Lackey and Wade [5^. The latter authors 
modeled the EOS as piecewise polytropes, allowing them 
to directly arrive at statements on the measurability of 
pressure as a function of density and neutron star radius 
as a function of mass. The latter method also has the ad¬ 
vantage that physical priors such as causality can more 
easily be folded in. Both [5S| and [25] concluded that 
A (mo), with mo = 1.4 Mq, could be measured with an 
accuracy of ~ 10 % by combining information from a few 
tens of sources. 

Of necessity, the studies in [531 US] used relatively 
simple waveform approximants, as otherwise the simu¬ 


lated data analysis problem would have been intractable 
with existing methods and computational infrastructure. 
Much effort is being put into large-scale numerical sim¬ 
ulations of the spacetimes of coalescing BNS, especially 
of the late inspiral P71l32j . The resulting waveforms are 
“hybridized” by matching them onto post-Newtonian or 
effective one-body waveforms, so that the earlier inspi¬ 
ral is also represented. While such waveforms represent 
the state of the art in our understanding of BNS coales¬ 
cence, producing a single one of them can take weeks. 
By contrast, high quality parameter estimation requires 
millions of waveforms to be compared with the data (see 
[33] and references therein). A full solution of the prob¬ 
lem of inferring the EOS from BNS detections will likely 
involve a combination of constructing phenomenological 
or “tuned” waveform models with input from numerical 
relativity [341140] , and significantly speeding up the anal¬ 
ysis of the data, e.g. through the use of Reduced Order 
Modeling; see [m HU and references therein. In that 
regard we note the recent work by Bernuzzi et al. [43] . 
who derived an effective one-body model that accurately 
describes tidal effects close to merger for a number of 
different EOSs, matching results from numerical simula¬ 
tions essentially to within the numerical uncertainties. 

On the other hand, when focusing on the inspiral 
regime, since the way tidal effects contribute to the phase 
is analogous for all of the PN approximants and the ef¬ 
fective one-body waveforms [numin], one might think 
that existing waveforms would already suffice to reliably 
extract information on the EOS from this part of the 
signal. However, as pointed out in [18] and studied in 
detail in [26l l46ll48j . significant biases can arise in the 
estimation of EOS effects due to discrepancies between 
waveform approximants - and presumably between these 
approximants and the true signal waveform - at high 
frequencies. Much of this is due to the fact that for 
the underlying point particle waveforms, the different PN 
waveforms and the effective one-body models differ sig¬ 
nificantly from each other at frequencies / > 400 Hz, 
where tidal effects become apparent. 

An important observation was made by Read et 
al. [3T], who studied the “distinguishability” ||(5Ii|| = 
a/(/ i 2 — hi\h 2 — hi) in terms of the usual PSD-weighted 
inner product (-I-) for waveforms hi, /12 of the same 
family but differing in their parameter values, in this 
case A. As can be seen in their Fig. 12, the dependence 
of ||(51i|| on changes in A is very similar for PN approx¬ 
imants and for hybridized numerical waveforms. Thus 
one may anticipate that PN approximants will allow us 
to predict how well one will be able to infer the EOS from 
GW measurements when sufficiently accurate waveforms 
will eventually become available for use in data analysis 
algorithms, even if the latter waveforms and appropri¬ 
ate analysis techniques are not yet at our disposal to¬ 
day. This will then inform the waveform modeling and 
data analysis communities as to what can reasonably be 
expected in terms of scientific output once their consid¬ 
erable efforts have come to fruition. Providing such an 
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assessment is the aim of the present paper. 

In this work we significantly expand on the results 
of [251126] in two ways: (a) we provide much more ex¬ 
tensive statistics through analysis of a larger number of 
simulated BNS sources, so as to identify worst-case and 
best-case scenarios depending on the detector noise re¬ 
alizations and the astrophysical distribution of the de¬ 
tected sources; and (b) we take into account as many 
physical effects as have been modeled, including tidal ef¬ 
fects to highest known order US], neutron star spins, the 
quadrupole-monopole interaction [491 ED] , the impact of 
early waveform termination due to the finite neutron star 
radii, and a relatively strongly peaked Gaussian distri¬ 
bution for the neutron star masses that is expected on 
astrophysical grounds [5X1454] . The effect of the latter 
was ignored in previous work, but as we shall see, it has 
a significant impact on inference of the EOS. As in |5S], 
we use two Bayesian data analysis methods: hypothesis 
ranking within a set of possible EOSs to see which one the 
true EOS is closest to, and parameter estimation on co¬ 
efficients in a Taylor expansion of the tidal deformability. 
We find that when source component masses are in a nar¬ 
row distribution, using a flat component mass prior for 
the analyses (as in previous work) increases the number 
of detections needed to distinguish between a soft, mod¬ 
erate, and stiff equation of state, and causes additional 
biases in parameter estimation on top of the ones due to 
imperfect knowledge of the signal model. On the other 
hand, this would get mitigated if we could assume accu¬ 
rate knowledge of the astrophysical mass distribution for 
neutron stars in binaries, so that it could be used as the 
prior distribution of masses. 

This rest of this article is structured as follows. In 
Sec. nil we introduce the waveform model and the EOS- 
related contributions from the effects mentioned above. 
In Sec. in we explain the two main methods used in the 
simulated data analysis: hypothesis ranking and param¬ 
eter estimation. Sec. jlVj explains the setup of our simu¬ 
lations. In Sec.jVjwe show the main results of this paper. 
A summary and discussion is given in Sec. m Finally, 
in the Appendix we further investigate the impact of the 
prior on component masses. 

Throughout this paper we will use units such that G = 
c = I unless stated otherwise. 


II. WAVEFORM MODEL AND EFFECTS OF 
THE NEUTRON STAR EQUATION OF STATE 

In this section we first discuss the general form of our 
waveform model, and then the way in which EOS effects 
enter. 


approximation (SPA), which yields a convenient analytic 
expression of the observed GW strain in the frequency 
domain [551 [5S] : 


AI, ?7) j2/3 

( 1 ) 

Here D is the distance to the source; {6, (p) denote the 
sky position with respect to the interferometer; (i, ip) de¬ 
termine the orientation of the orbital plane in relation 
to the observer; M = is the chirp mass, with 

M = mi -I- m 2 being the total mass and rj = mim 2 /M^ 
the symmetric mass ratio; tc and are, respectively, 
the time and the phase at coalescence; and xi, X 2 are 
the neutron stars’ dimensionless spins. The “frequency 
sweep” F{f-M,r],xi,X 2 ) (related to the time domain 
phase $(t) by F{f) = $(t(/))/7r) is an expansion in 
powers of / with coefficients that depend on masses and 
spins, and the SPA phase takes the general form 


4'(/) = 27r/4 -pc - x + 
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where the tpk and ^ again depend on masses and spins. 
For the low-mass systems considered in this paper, the 
Advanced LIGO-Virgo network will not be very sensitive 
to subdominant PN contributions to the amplitude [ 571 - 
[59] . and we use the “restricted” post-Newtonian approx¬ 
imation, in which only leading-order PN contributions to 
the amplitude are taken into account. In particular, this 
means that spin and EOS effects appear in the phasing 
only. 

For the purposes of this paper, the phase was taken to 
3.5PN order, with inclusion of spin effects up to 2.5PN 
following ED], and we refer to that paper for explicit ex¬ 
pressions. For simplicity we assume that the components’ 
spins are aligned or anti-aligned with each other and with 
the direction of orbital angular momentum; at the time 
this work was started, frequency-domain, precessing-spin 
waveform approximants like the ones of Lundgren and 
O’Shaughnessy EH, of Hannam et al. El iO], and of 
Klein et al. [62] were not yet available. It is quite pos¬ 
sible that inclusion of precession would aid in breaking 
degeneracies between spins and mass ratio, as was sug¬ 
gested in e.g. E3]j enabling more accurate measurements 
of EOS effects. 

We take into account three ways in which the EOS af¬ 
fects the waveform: tidal deformations, the quadrupole- 
monopole effect, and the possible early termination of 
the waveform due to the finite size of the neutron stars, 
whose radii are set by their masses and the EOS. Let us 
discuss these in turn. 


A. General form of the waveform model 


B. Tidal deformations 


We model gravitational waveforms from the quasicir¬ 
cular inspiral of BNS systems using the stationary phase 


Towards the end of the evolution of a BNS system, 
when the gravitational wave frequency reaches / > 400 
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Hz [T5], the tidal tensors £ij of one component of the 
binary will start to induce a significant quadrupole mo¬ 
ment Qij in the other. In the adiabatic approximation, 
the two are related by [HJ |M1 [SS] 

Qij = (3) 

where m is the mass of the neutron star that is experienc¬ 
ing the quadrupole deformation, and the function A(m) is 
the tidal deformability, which is determined by the EOS. 
The deformations of the two neutron stars in turn affect 
the orbital motion, and this is one way in which the EOS 
gets imprinted upon the gravitational waveform. The 
deformability A(m) is related to the second Love number 
k 2 (jn) and the neutron star radius R{m) through A(m) = 
(2/3) fc 2 (m) Tidal effects only enter the phase 

starting at 5PN order [5S], but as mentioned before, the 


prefactors are sizable (A/M® oc (i?/M)® ~ 10^ — 10®), 
which is why we can hope to infer information on the 
EOS from the tidal deformation. 

The effects of tidal deformations on the orbital motion 
were calculated up to IPN (or 6 PN in phase) by Afines, 
Flanagan, and Hinderer [li], and more recently to 2.5PN 
(or 7.5PN in phase) by Damour, Nagar, and Villain |19) . 
The latter expression is what we will be using in this 
paper; for completeness we reproduce it here. In terms 
of the characteristic velocity v = (ttM/)^/®, one has 

^l[v) = Tpp(u) -I- 4'tidal(?^), (4) 

where 4'pp(p) is the phase for the inspiral of point parti¬ 
cles, and 4'tidai(^^) is the contribution from tidal effects. 
The latter takes the form 


4'tidai('y) = 
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where Xa = tua/M, A = 1,2, and Aa = X{mA)- We 
should note that the calculation leading to this expression 
ignores (i) contributions from higher-order multipoles as 
these are estimated to give small corrections, and (ii) a 
number of as yet unknown functions that appear in the 
7PN phase contribution; in m these too were argued 
to be negligible and we refer to that paper for details. 
Contributions to the phase at increasing PN order, for a 
BNS system of (1.35, 1 . 35 )M 0 with a stiff (MSI) EOS, 
are illustrated in Fig. 4. 

For the function X(rn), in our simulated signals we will 
use quartic polynomial fits to predictions corresponding 
to different EOS from Hinderer et al. |18| . with maximum 
residuals of ~ 0.02 (which will turn out to be negligible 
compared to the measurability of A). Examples of such 
fits for a soft (labeled SQM3), a moderate (H4), and a 
stiff EOS (MSI) are shown in Fig. 1. 


C. Quadrupole-monopole effects 

As mentioned before, tidal effects are not the only way 
the EOS enters into the gravitational waveform. If a 
neutron star is spinning, it takes on an oblate shape. As¬ 
suming an axisymmetric mass distribution with respect 
to the axis of rotation, the deformation can be expressed 


to leading order by means of a dimensionless quadrupole 
moment parameter g, defined as |50] 

J *"(A6')^2(cos6»)dcos6i, (6) 

where P 2 {x) = {3x^ — l)/2 is the second Legendre poly¬ 
nomial, and 1 / is a potential related to the metric of 
a stationary axially symmetric body; more specifically, 
the line element in the form introduced by Komatsu- 
Eriguchi-Hachisu [BB] reads: 

ds^ = -I- sin^ 9 {d4> - ojdtf 

+ 6 ^“ (dr^ -h r^dO"^) , (7) 

where the undetermined a, j3, v are all functions of (r, d). 
The quadrupole moment q is the leading-order (1/r®) 
coefficient of the second multipole in the asymptotic ex¬ 
pansion of v{r, 9) and can be calculated numerically. This 
quantity is the general-relativistic equivalent of the New¬ 
tonian mass quadrupole moment. 

Since a stiffer EOS implies a larger neutron star (NS) 
radius for a given mass, the quadrupole moment increases 
in absolute value with the stiffness of the EOS. Examples 
of q estimates for different EOS were calculated numer¬ 
ically in m based on the expressions of Ryan [smsH]. 
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FIG. 1: The tidal deformability parameter A(m) as a func¬ 
tion of neutron star mass for three different EOS: a soft 
one (SQM3), a moderate one (H4), and a stiff one (MSI). 
Adapted from |18| . Curves are fitted quartic polynomials, 
whose residuals are shown in the lower subplot. Only masses 
within the unshaded region [1, 2]Mq will be considered in our 
analyses. 



FIG. 2: The quadrupole parameter a(m) as a function of neu¬ 
tron star mass for the three different EOSs in Fig. 1. The hor¬ 
izontal dashed line indicates the value for black holes, which is 
a = 1 m- Only masses within the unshaded region [1, 2\Mq 
will be considered in our analyses. 


depends on masses and spins through 


These demonstrated the dependence on the dimension¬ 
less spin X, which for a fixed NS mass can be fit very well 
up to the maximum spin value Xmax (also dependent on 
the EOS) by a quadratic rule: 


CTQM 


- H 


QA 


A=l,2 


/ruAy 
\ M J 


3(xa - Lf-l 


( 10 ) 


«(^t) (^) 3(xa • L )^ 


- 1 


A=l,2 




q ~ -ax^ (8) 

where a = aEOs(tn) is a mass-dependent parameter. Fur¬ 
ther evidence to support the quadratic relation Eq. § is 
given in [Si [70]. The authors of [Sim] also point out a 
spin correction in the identification of multipole moments 
that was previously overlooked; this correction preserves 
the quadratic spin behaviour of Eq. ([^, and vanishes in 
the slow-rotation limit. Assuming that this relation will 
hold for any EOS, we will only be concerned with the 
spin-independent parameter a which, similar to the tidal 
deformability parameter A, has a functional dependence 
on the neutron mass that is determined by the EOS. 

The effect of such a quadrupole moment on the grav¬ 
itational waveform emitted by a binary system was de¬ 
rived in |49| . To Newtonian order it introduces an ad¬ 
ditional coupling in the effective gravitational potential, 
between the mass quadrupole of each spinning neutron 
star and the mass of its companion, whence the name 
“quadrupole-monopole (QM) effect”. In the stationary 
phase approximation, the additional contribution to the 
GW phase due to the QM interaction reads: 

T / ^ 30 

^'QM(^^) =, (9) 

making it of 2PN order in phase. The parameter ctqm 


where the unit vectors xa are the direction of the spins. 
In the last line we used the rule (|^; we see that with 
this assumption, 'I'qm(^^) is quadratic in the component 
spins. Finally, note that in the case of (anti-)aligned 
spins, which we will assume throughout, 3(xa • — 1 = 

2 . 

As mentioned above, in our simulations we will use 
predictions for A(m) corresponding to different EOSs 
from m- In order to compute a{m), we make use of 
the recently discovered phenomenological Love-Q rela¬ 
tion [zuizg, which is believed to hold irrespective of the 
EOS: 


In a(TO) = 0.194-b 0.0936 In ^-f 0.0474 (In 




m" 


-4.21 X 10"^ ( In^^ -b 1.23 X lO"'^ fin ^ 

' -.b ] \ 
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( 11 ) 


The relative fractional errors due to the universal fit were 
estimated in m for several EOSs to be at the 1% level. 
Together with Eqs. (101 and (|^, this then allows us to 
compute the QM contribution to the phase. Fig. 2 shows 
a{m) for the EOSs in Fig. 1. QM contributions to the 
phase are expected to be subdominant compared to the 
tidal effects of Sec. |IIB[ even for relatively fast spinning 
NS, as shown in Fig. 4. 
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FIG. 3: The frequencies /lso and /contact as functions of mi, 
m 2 for the EOS shown in Fig. 1. 



f[Hz] 

FIG. 4: Phase contributions of the QM effect and tidal effects 
up to different PN orders as functions of GW frequency for 
a (1.35,1.35)M0 binary with a stiff EOS (MSI). The QM 
contribution from each NS scales quadratically with its spin 
and is shown here for xi = Xa = 0.1. The dashed vertical 
lines indicate the contact and LSO frequencies. 


D. Termination of the waveform at contact 


In the recent simulations |25l I26j . the waveform was 
cut off at a frequency corresponding to the last stable 
circular orbit (LSO) in the point particle limit, given by 


/lso = 


1 


6^/^7rM 


( 12 ) 


However, as we shall see below, it will often happen 
that the two neutron stars attain physical contact be¬ 
fore the corresponding distance between the components 
is reached. In this paper, we instead impose the cutoff 


/cut — Cninl /lso )/contact} 5 


(13) 


where, using Kepler’s third law, the “contact frequency” 
is given by 


/< 


contact — 


M 


TT \R{mi) + R{m 2 ) 


1/2 


(14) 


We stress that the termination condition © is still 
a heuristic one, but it will be more realistic than termi¬ 
nation at /lso- Moreover, the length of the waveform 
itself carries physical information m, in this case on the 
EOS, which we wish to incorporate [M]. On the other 
hand, shorter waveforms have a smaller number of cycles 
from which information can be extracted; when we come 
to the results of our simulations we will see which effect 
wins out. 

In order to compute the radii R{mi), R{m 2 ), we again 
make use of a recently discovered phenomenological re¬ 
lation, this time between the compactness C = m/R and 
A |7B]: 


C = 0.371 - 3.91 X 10"^ In A -f 1.056 x 10 


-3 


In 


A 


m" 


( 15 ) 


For a given EOS (ie a given relationship A(to)), the 
above expression gives us R{m), from which the con¬ 
tact frequency (141 is obtained. The relative error in 
the compactness (and hence in the radius) due to the fit 
of Eq. (15) was found to be at the 2% level, implying a 
similar error in the contact frequency. 

Fig. 3 shows the dependence of /lso and /contact on 
component masses mi, m 2 for the EOS considered above. 
Note how in the astrophysically relevant range mA S 
[1, 2 ] Mq, H = 1 , 2, it often happens that /contact < /lso, 
especially for low masses and for the stiffer EOS (MSI) 
which can support larger neutron star radii. 


III. BAYESIAN METHODS FOR INFERRING 
THE NEUTRON STAR EQUATION OF STATE 


In this section we present two qualitatively different 
Bayesian methods that one may use to acquire informa¬ 
tion on the neutron star equation of state: (i) hypothesis 
ranking for different proposed EOS based on how well 
each of them matches the available data, and (ii) the es¬ 
timation of parameters which for a given EOS will be the 
same across sources. Both of these allow us to combine 
information from multiple detections so as to arrive at a 
stronger result. These methods were already explained 
in [5S]; for completeness we recall the basic ideas. 


A. Hypothesis ranking 

Given a set of (finitely many) EOS models 
{Ml, M 2 ,..., Mk}, we will be interested in ranking them 
in the light of the available data. The ranking process 
will be on a set of hypotheses {Hi-,i = ,K{, where 

Hi states that Mi is the true model for the neutron star 
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EOS. Each model comes with a particular functional de¬ 
pendence of the tidal deformability on the mass, 
and hence a QM parameter an d ra dius- mas s re¬ 
lation obtained through Eqs. (Ill and (15), re¬ 

spectively. The general form of the waveform given by 
Eqs. ([^ and ® , together with the tidal and QM contri¬ 
butions ([^, ([^ and the waveform termination condition 
(13), yields a waveform model h^^\f;6) associated with 


the hypothesis Hi. The parameter spaces {9} (masses, 
spins, sky position, orientation, distance, time at coa¬ 
lescence, and phase at coalescence) are the same for all 
the hypotheses Hi, i = but for given com¬ 

ponent masses itia, A = 1,2, the calculated values for 
Ay 4 = A^^^ttia), CIA = a^''\mA), and differ between 
one hypothesis and the next. 

Given a detection d (to be thought of as a stretch of 
detector data containing a confirmed BNS signal), a hy¬ 
pothesis Hi, and any background information / that we 
may have, the likelihood function is defined as [TTHzg 


p{d\H„e,I) 

^ SM) 


= Af exp 


, ( 16 ) 


with Sn{f) being the detector’s one-sided noise power 
spectral density (SO]; d{f ) is the Fourier transform of the 
data stream, and A/” is a normalization factor. We will 
take the low-frequency cutoff to be flow = 40 Hz, and 
/cut is the one in Eq. (13). To compute p{d\Hi,9,I) we 


will use the method of nested sampling as implemented 
in GW parameter estimation by Veitch and Vecchio [771- 
ES]; see also [33]. The evidence is given by 


P{d\H,,I)= / dep{9\I)p{d\H,,9,I), 


(17) 


with p{9\I) being the prior density distribution. Using 
Bayes’s theorem, the posterior probability of the hypoth¬ 
esis Hi given the data d is then obtained through 


P{H,\d,I) 


P{d\H„I)P{H,\I) 

P{d\I) 


(18) 


where P{Hi\I) is the prior probability for Hi before any 
measurement has taken place, and P{d\I) is the prior 
probability of the data. Finally, the odds ratio between 
any two hypotheses Hi, Hj is defined as 


, ^ Pmd,I) ^ P{H,\T) P{d\H,,I) 

^ P{H,\d,I) P{H,\I)P{d\H„iy 


Note that the unknown P{d\I) conveniently drops out of 
this expression. 

The above framework can trivially be generalized to 
the case of multiple detections di,d 2 , ■ ■ ■ ,dN■ Using the 
multiplication rule for independent random variables, one 
obtains 


{N)^^ ^ Pmi) A Pidn\H,^I) 

^ p(i7,i/) n p(d„|i7„/)- 


( 20 ) 


The probabilities P{Hi\I) quantify our prior degree of 
belief in the hypotheses. Gurrently a wide range of 
EOS are still consistent with existing observations m, 
including the ones that we will use in our simulations 
(although general theoretical considerations suggest a 
more restricted range; see e.g. [ST]). In the absence of 
any additional information, it then makes sense to set 
P{Hi\I) = P{Hj\I) for all i and j. 

If > 1 , then the data favors the hypothesis Hi 

over the hypothesis Hj . By looking at the odds ratios for 
all pairs of hypotheses, we arrive at an overall ranking of 
all the Hi,H 2 , ■ ■ ■, Hk- We explicitly note that 

1. Even if the true equation of state were in the set 
Hi, i = 1, ..., K , one should not necessarily expect 
it to end up at the top of the ranking; this is due to 
the effects of noise and the fact that the majority of 
detected sources will have low signal-to-noise ratios 
(SNRs). 

2. In practice, the correct equation of state will prob¬ 
ably not be in the finite set Hi, i = 1,..., K . Nev¬ 
ertheless, one may expect the highest-ranked hy¬ 
pothesis to be close to the true one. 

Here a notion of closeness or distance in a space of func¬ 
tions is implied; this can be defined by e.g. employing the 
norm ||/|| = (J |/pd^)^/^. The integration measure p 
need not be uniform in mass {i.e. in principle dp y dm), 
but should rather reflect the amount of information that 
is collected from each infinitesimal mass interval. That 
is, if two functions differ significantly at a mass inter¬ 
val where no sources are found, but are almost equal 
elsewhere, then the “distance” between them should be 
small. Here however, the set of functions A*-®^(to) that 
we consider are clearly distinguishable across the mass 
interval of interest [ 1 , 2 ]M 0 and admit a strict ordering 
in terms of stiffness. 

Finally, we note that it is often convenient to work 
with the logarithms of the odds ratios, as we will also do 
here. 


B. Parameter estimation 


An obvious advantage of hypothesis ranking is that in¬ 
formation from multiple detections can trivially be com¬ 
bined; see Eq. (20). In measuring parameters, we will 


want to do the same. Hence it will not do to just 
measure e.g. the tidal deformabilities Aa, A = 1 , 2 for 
each source individually, since these numbers are mass- 
dependent and will vary from source to source. In order 
to combine information across sources, one must identify 
observables that only depend on the equation of state 
and not on incidental details of the sources. In [25], the 
function A(m) was approximated by a Taylor expansion 
around some reference mass mo, truncated at a suitably 
















high order: 


Jmax -1 / 

-^1 / TO - Too 

Mm) = > — c,- ——- 

V ^0 


( 21 ) 


For a given EOS, the coefficients cj will be the same for 
all the sources. 

A waveform model h{f;9,{cj}) can be constructed 
along the lines of Sec. |l^ but this time substituting the 
expansion (21) for A(to); here 9 still represents masses, 
spins, sky position, orientation, distance, time of coales¬ 
cence, and phase at coalescence. Given a detection d, a 
posterior density function for each of the Cj can be con¬ 
structed [25l [77U79] . For instance, in the case of cq. 


p{co\d,I) = / d 6 »dci...dcj„^^p( 6 », {cjlld,/), (22) 


where the joint posterior density function for all the pa¬ 
rameters takes the form 


p{9,{c,}\d,I) 


p{d\9,{cj),I)p{9,{cj)\I) 

p{d\I) 


(23) 


Here l/p{d\I) acts as a normalization factor, and 
p{9, {cj}\I) is the joint prior density of all the param¬ 
eters; we will assume that the latter can be written as 


p{9,{c,}\I)=pm np(c,|/), (24) 

j=o 

where p{9\I) and p{cj\I) are separate prior densities for 
the 9 and all the Cj, respectively. Finally, the likelihood 
p{d\9,{cj}J) is given by 


p{d\9,{cj},I) =J\f exp 



d/ 


\dif)-hif-9,{cj})\^ 

Snif) 


(25) 


Given multiple detections di, c? 2 ,..., dAr, individual 
posterior density functions such as p(co|d„,/) can triv¬ 
ially be combined [ 25 ] : 


N 

p(co|di,d2,...,dAr) =p{co\I)^~^ 

n—1 


where we again assumed independence of the d„ and have 
used Bayes’s theorem. Similar expressions can of course 
be obtained for p(cj |di,d 2 ,... ,djv), j = 1 , ■. ■, Jmax- 
A choice needs to be made for the order jmax at which 
the Taylor expansion (21) is truncated. In [25], the au¬ 
thors chose jmax = Ij Under the expectation that but 
two coefficients will be measurable when EOS effects only 
enter the waveforms through the two parameters A(toi), 
A(to 2 ). Here we will instead use a quadratic approxima¬ 
tion to A (to): 


A(to) ~ Co -I- Cl 


TO — Too 

Mr. 


1 

+ 2C2 


m — mo 

Mr. 


(27) 


with Too = 1.4 Mq. Visual inspection of the A(to) for 
the EOSs considered in e.g. [H] (see their Fig. 2) already 
suggests that in the most plausible mass range for neu¬ 
tron stars (roughly to G [ 1 , 2 ]M 0 ), this will tend to be 
a good approximation, which is why we make the choice 
here; see also Fig. 1 for the EOSs used in this paper. As 
we shall see, in practice neither ci nor C 2 will be mea¬ 
surable even with a large number of sources, but cq will 
be. 


IV. SETUP OF THE SIMULATIONS 

We now briefly describe how the simulations were set 
up. Different choices were made for the parameter dis¬ 
tribution of the simulated signals, or injections, and for 
the priors on the parameters used for the data analysis. 


A. Injections 

For the parameters of the simulated sources, we choose 
astrophysically motivated distributions. Sources are 
placed uniformly in comoving volume in a distance range 
D G [100, 250] Mpc. The upper bound is approximately 
the angle-averaged range which Advanced LIGO is ex¬ 
pected to reach, while the lower bound corresponds to 
the distance within which one can expect a detection 
once every two years [7]; this cutoff also serves to ex¬ 
clude unexpectedly loud sources from our ensembles. The 
sky location (9, (p) and orientation (t, tf) are both dis¬ 
tributed uniformly on the sphere. The phase at coa¬ 
lescence ifc is taken to be uniform in [0, 27r). For the 
spins we note that observed pulsar periods and assump¬ 
tions about spin-down rates lead to birth periods in the 
range 10 — 140 ms [82] , which corresponds to dimension¬ 
less spins J/to^ < 0.04; the fastest known pulsar in a 
BNS system has J/to^ ~ 0.02. The observed popula¬ 
tion is expected to spin down to much lower values of 
spin at the time of coalescence, but currently there is no 
good estimate on the spin distribution for a population 
of coalescences that will be observed with GW detectors. 
We choose to take a conservative approach and assume 
the spins to be small but non-negligible for the analysis: 
when including spins in the simulated signals, we take 
them to be Gaussian distributed with zero mean and a 
spread of = 0.02. Unless stated otherwise, the com¬ 
ponent masses are picked from a Gaussian with mean 
firn = 1.35Mq and spread = 0.05 Mq. The latter is 
inspired by estimates of the mass distribution of known 
neutron stars in BNS systems [5TI - I55] . 

For the EOS, we want to find out under what circum¬ 
stances one will at least be able to distinguish between 
a stiff, a moderate, and a soft EOS. For these we choose 
the EOSs labeled MSI, H4, and SQM3, respectively, in 
Fig. 1. 

The simulated GW waveforms are added coherently 
to simulated data streams for Advanced LIGO detectors 












9 


at Hanford, WA, and Livingston, LA, as well as an Ad¬ 
vanced Virgo detector at Cascina, Italy. The detector 
noise is taken to be stationary and Gaussian, where the 
underlying noise curve for Advanced LIGO is the one 
with zero detuning of the mirror and high laser power, 
and for Advanced Virgo we choose a BNS-optimized 
curve assuming an appropriate choice for the signal re¬ 
cycling detuning and the signal recycling mirror trans¬ 
mittance; see [5S] and references therein. To realisti¬ 
cally simulate a scenario in which BNS signals have been 
detected, we impose two conditions on the signals: (i) 
the optimal network SNR should be greater than 8 but 
smaller than 30, and (ii) the post-analysis log Bayes fac¬ 
tor for signal versus noise should be greater than 32 with 
templates that assume point particle coalescence [55] . 
This is also what was done in and we refer to that 
paper for more details. 


sider MSI, H4, and SQM3, as well as the “point particle” 
model, denoted PP. The latter corresponds to a waveform 
model where A(m) = 0. 

Finally, when doing parameter estimation on the co¬ 
efficients Co, Cl, and C 2 in the quadratic approxima¬ 
tion to A(m) as in Eq. (27), the priors are chosen to 
be Co e [0,5] X 10-23s^ ci e [-2.5,0] x lO-^SgS, and 
C 2 G [—3.7,0] X jn the mass regime of interest, 

this captures all the EOSs in Fig. 2 of m 


V. RESULTS 


Let us now present the results of our simulations, first 
for the hypothesis ranking described in Sec. Ill A and 


then for parameter estimation as explained in Sec. IIIB 


B. Templates 

The data analysis was performed as described in 
Sec. |III[ with the following choices of prior distributions 
for the parameters. Distance was allowed to vary in the 
range D S [1,1000] Mpc. The coalescence phase (pc was 
taken to be uniform on [0,27r). The coalescence time was 
allowed to vary within 100 ms of an injection. When al¬ 
lowing for non-zero spins in the templates, we took the 
priors on xi, X 2 to be uniform on the interval [—0.1, 0.1]. 

As we shall see, the prior density distribution for the 
component masses will play an important role. In princi¬ 
ple we could take this to be the same as the mass distribu¬ 
tion for the injections, i.e. Gaussian with /r = 1.35 Mq 
and dm = 0.05 Mq. However, we would then implic¬ 
itly be assuming that the astrophysical mass distribu¬ 
tion of neutron stars in binaries will be reliably known 
in the advanced detector era. At the time of writ¬ 
ing only 9 double neutron star systems have been ob¬ 
served, sometimes with large error bars on the mea¬ 
sured masses; it seems unlikely that this situation will 
improve dramatically in the next few years. We also 
note the differing results for observationally based esti¬ 
mates of the mass distribution in BNS systems; for ex¬ 
ample, (p-mjCTm) = (1.37 Mq, 0.042 Mq) in Valentim et 
al. El] and {p.m,crm) = {1.33 Mq, 0.13 Mq) in Kiziltan 
et al. [ 55 , the difference partially being due to the use 
of different subsets of the known systems based on the 
reliability of individual mass measurements. Finally, it 
is possible that due to selection biases, the distribution 
of masses in electromagnetically observed neutron star 
binaries will not be identical to the mass distribution 
in BNS coalescences detected by Advanced LIGO and 
Virgo. For these reasons, we will mostly assume a flat 
component mass prior with m £ [1,2]M0. However, in 
the Appendix we will also briefly investigate what hap¬ 
pens if the astrophysical distribution of masses of neutron 
stars in binaries can be assumed known after all. 

For the EOSs in the hypothesis ranking, we again con- 


A. Hypothesis ranking 


A first estimate of how well one will be able to deter¬ 
mine the EOS using hypothesis ranking was presented 
in our earlier paper |25j . In that work, only tidal ef¬ 
fects up to IPN rather than 2.5PN were taken into ac¬ 
count, quadrupole-monopole contributions were disre¬ 
garded, waveforms were terminated at LSO instead of 
the frequency /cut of Eq. (13), and spins were set to zero 
both in the injections and in the template waveforms. 
Additionally, the component masses were taken to be dis¬ 
tributed uniformly in the interval [1,2]Mq rather than 
according to a Gaussian with mean fj,m = 1.35 Mq and 
spread = 0.05 Mq. 

Ideally one would like to look at the impact of each 
of these effects individually. However, the simulations 
presented in this paper are computationally expensive if 
one wants to have good statistics. For this reason, we 
proceed as follows: 


• First we set the spins to zero both in injections and 
templates (so that the quadrupole-monopole effect 
is not present), but we take tidal effects to 2.5PN 
and terminate the waveforms at the minimum of 
contact frequency and LSO frequency. We generate 
results for injected component masses distributed 
uniformly in [1,2]Mq, and then for component 
masses following a Gaussian with = 1.35 Mq 
and dm = 0.05 Mq; however, in both cases the 
mass prior in our Bayesian analysis is taken to be 
uniform on [1,2]Mq. Again because of computa¬ 
tional cost, we only make this comparison for the 
case where the EOS in the signals is MSI, i.e. the 
stiffest equation of state considered in this paper. 


• Next we specialize to the more astrophysically mo¬ 
tivated Gaussian distribution for the component 
masses (still keeping a uniform prior in the analy¬ 
sis), and we also switch on spins. In the injections, 
we let the latter be Gaussian distributed with zero 
mean and d^ = 0.02, while in the templates we let 











10 




FIG. 5: Hypothesis ranking in the case where the EOS of the simulated sources is MSI. As in [^, spins are set to zero both 
in the injections and the templates, and masses are distribnted uniformly on the interval [1,2] Mq, so that any qualitative 
differences with previous results come from considering tidal effects to higher PN order, and terminating waveforms at contact 
or LSO, whichever comes sooner. Left: The cumulative distributions of the log odds ratios In for catalogs of 20 sources 

each, where “EOS” is in turn H4, SQM3, and PP. Note how the EOSs are ranked according to how dissimilar they are to the 
correct one: PP differs the most and is indeed the most deprecated, followed by SQM3 and H4. Right: The fraction of catalogs 
for which PP, SQM3, and H4 are correctly ranked lower than MSI, as a function of the number of events per catalog. What is 
shown are the medians and 95% confidence intervals obtained from combining individual sources into catalogs in 1,000 different 
ways. 



InOj^sf # events / catalogue 


FIG. 6: The same as in Fig. 5, but this time with a relatively strongly peaked Gaussian distribution for the injected component 
masses (while still using a flat mass prior in the analyses). To approach the discemibility of EOS seen in Fig. 5, we now need 
0(100) sources per catalog. Even then, H4, the EOS that is closest to the injected MSI, can not be distinguished from it. 


the prior on the spins be uniform on the interval 
[—0.1,0.!], to reflect the ignorance about spins we 
will in practice have. Since in this case we are in¬ 
cluding all the astrophysical effects considered in 
this paper, we generate results not only for MSI, 
but for H4 and SQM3 as well. 


1. Zero spins; flat versus Gaussian distribution of 
component masses 

First we consider the case of zero spins in injections 
and templates, and component masses are distributed 


uniformly on the interval [1, 2] Mq. Results are shown in 
Fig. 5. We let the injections have MSI as their EOS, and 
as in , we compute the log odds ratios In Omsi 
catalogs of 20 sources each, where, in turn, “EOS” stands 
for PP, SQM3, and H4. Examples of the cumulative dis¬ 
tributions of these log odds ratios are shown in the left 
panel of the figure. In the absence of detector noise, one 
would have In^^OoEOS < q in all three cases, since any 
EOS different from the correct one (MSI) would be dep¬ 
recated. What we see is that In^^^^O^g^ < 0 for about 
80% of the catalogs, while < 0 in about 60% 

of the cases. Note that H4 is the most similar to MSI, 
followed by SQM3 and PP; and indeed, the log odds ra- 
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tios obtained tend to correctly rank the EOSs in this 
way. This is similar to what one sees in the top right 
panel of Fig. 2 of However, despite the fact that 
in the present work we take tidal effects to much higher 
order, the left tails of the cumulative log odds ratio dis¬ 
tributions stretch to less negative values. This can be 
explained by the different termination of the waveforms, 
which for the EOSs and mass distributions we consider 
tends to be at contact rather than LSO (see Fig. 3). For 
a typical system with component masses (1.35,1.35) Mq 
and equation of state MSI, the termination frequency is 
/contact = 1222 Hz whereas /lso = 1629 Hz, so that the 
signal contains less information on tidal effects (which 
manifest themselves at high frequency) than in [25] . In¬ 
deed, as can be seen in Fig. 4, the higher-order tidal 
effects due to their alternating signs do not significantly 
change the number of cycles in the phase, though they 
will add some structure because they come with different 
powers of u; on the other hand, termination at contact 
seems to have a much stronger effect, cutting the tidal 
phase short (in this example by roughly 15 rad). The 
QM effect is much weaker and is not expected to give a 
significant contribution to the inference. 

The right panel of Fig. 5 shows the fraction e^Qg of cat¬ 
alogs for which MSI is ranked higher than, respectively, 
H4, SQM3 and PP (f.e. In^^ofoEOS ^ q where “EOS” is, 
in turn, PP, SQM3, and H4), as a function of the number 
of sources per catalog. Also shown are 95% intervals on 
^EOS obtained from combining 1,800 individual sources 
into catalogs in 1,000 different ways. We see the same 
trend as in the left panel: H4, being the most similar to 
MSI, is ranked below MSI the least often, and PP, be¬ 
ing the most dissimilar, the most often. We note that in 
going to a higher number of sources per catalog, we start 
experiencing small number statistics; at 100 sources per 
catalog only 18 independent catalogs can be composed. 

Next, in Fig. 6 we look at the case where the spins 
are still zero in injections and templates, but the in¬ 
jected masses are distributed according to a Gaussian 
with = 1.35 Mq and am = 0.05 Mq. Unlike in Fig. 5, 
in the left panel showing the cumulative distributions 
of the log odds ratios, we now consider catalogs of 100 
sources each, which turns out to be necessary to approach 
the discriminatory power we had with a uniform mass 
distribution. Even then, H4, the EOS that most closely 
resembles the injected MSI, is not distinguishable from 
it: the probability that MSI gets ranked above H4 is ap¬ 
proximately the same as the probability that H4 ends up 
above MSI. 


2. Gaussian mass distribution, non-zero spins 

We now specialize to the astrophysically better moti¬ 
vated Gaussian distribution for the injected component 
masses (but sticking to a uniform mass prior in our anal¬ 
yses), and we switch on spins xAi A =1,2. In the injec¬ 
tions, the spins are Gaussian distributed with zero mean 


and cr^ = 0 . 02 , while in the templates, the XA have priors 
that are uniform on the interval [—0.1, 0.1]. This time we 
give results for injections where the EOS is MSI, H4, and 
SQM3, respectively. 

In the left panels of Fig. 7 we see examples of cumu¬ 
lative distributions of In for catalogs of 100 

sources each, where “inj” is the injected equation of state, 
while “EOS” is, in turn, taken to be each of the other 
three EOSs considered in this paper. From top to bot¬ 
tom, the injections follow MSI, H4, and SQM3, respec¬ 
tively. 

In the right panels of Fig. 7 we again vary the number 
of sources per catalog, and show the fraction 6^03 of 
times that the injected equation of state is ranked higher 
than each of the other three EOSs in turn. For a given 
number of sources per catalog, we combine individual 
sources into catalogs in many different ways and look at 
the medians and 95% confidence intervals of the egos- 

Let us first compare the results for MSI (top panels in 
Fig. 7) with the ones for Gaussian distributed masses but 
zero spins in injections and templates (Fig. 6 ). Looking 
at the e^os, we infer that EOSs again tend to be ranked 
correctly according to “stiffness” and similarity to MSI, 
and we even see some improvement in the discernibility 
of H4 from MSI, especially as the number of sources per 
catalog goes to 100 . 

For H4 injections (middle panels in Fig. 7), the me¬ 
dians of ejlos ordered, with the median of Cpp 

staying above that of Cmsi’ 'which in turn trumps e^qj^a- 
However, H4 being in between MSI and SQM3 in stiff¬ 
ness (see Fig. 1), the 95% uncertainty intervals of the 
^EOS show considerable overlap; although H4 is ranked 
above each of the other EOSs reasonably frequently, the 
internal ranking is less clear. 

Finally, for SQM3 (bottom panels), this being the soft¬ 
est EOS other than the PP model, the stiff MSI tends 
to be deprecated reasonably strongly, but it is hard to 
distinguish SQM3 from either H4 or PP. 


B. Parameter estimation 


We now turn to the data analysis setup described in 


Sec. HI B Here the templates used for the recovery do not 
have a fixed tidal deformability function A(to); rather, it 
is modeled by a quadratic polynomial as in Eq. (271, 


where the coefficients cq , ci, C2 are now free parameters 
to be estimated, on top of all the usual ones (masses, 
spins if applicable, time and phase at coalescence, sky 
position, orientation, and distance). To the extent that 
the quadratic approximation can capture the EOS in the 
signal in the relevant mass range, in the measurement 
process we can assume cq, ci, and C2 to have fixed values, 
so that their posterior densities can be combined across 
sources as in Eq. (26). 

In our earlier paper [25] . where only a linear approx¬ 
imation to A(m) was used, it was found that only the 
zeroth-order coefficient was measurable. The quadratic 





12 



InO 


EOS 

MSI 


# events / catalogue 





FIG. 7: The same as in Fig. 6, except that simulated sources have (anti-)aligned spins sampled from a Gaussian distribution 
centered at zero and with = 0.02. However, the prior on spins used in the recovery is uniform on the interval [—0.1,0.!]. 
Left panels: Examples of cumulative distributions of log odds ratios for the injected EOS versus the other ones considered, for 
catalogs of 100 sources each. From top to bottom the injected EOS is MSI, H4, and SQM3, respectively. Right panels: The 
fraction of catalogs for which the correct EOS is ranked higher than each of the others in turn, as a function of the catalog size. 
Here too, medians and 95% confidence intervals are shown, obtained from combining sources into catalogs in 1,000 different 
ways. 


approximation used in the present paper should allow 
for a better fit, but here too, it turns out that only the 
leading-order coefficient Cq can be measured with any 
kind of accuracy. Thus, unlike with hypothesis ranking, 
in practice only a single number pertaining to the EOS 
is being extracted from the data. Nevertheless, one has 


Co = A(mo), with mo some fixed reference mass (which 
we will take to be I. 4 M 0 ), and as can be seen in Fig. 2 
of the paper by Hinderer et al. |18j , which shows nearly 
20 different predictions for A(m), valuable information 
could be gleaned from just that one number. 

As before, we consider the following cases: 
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• Spins are zero both in injections and templates, 
and we compare results for an injected mass dis¬ 
tribution that is uniform on [1,2] Mq with what 
one gets with a Gaussian mass distribution that 
has Urn = 1.35Mq and am — 0.05M©. However, 
for the templates we do not assume knowledge of 
the astrophysical mass distribution, sticking to a 
uniform mass prior on [1, 2] Mq. 

• Next we specialize to the Gaussian injected mass 
distribution, and switch on spins. In the injected 
waveforms, the latter are drawn from Gaussian dis¬ 
tributions with zero mean and a^ = 0.02, while in 
the templates the priors for the spins are uniform 
on [—0.1,0.!]. 

We stress again that for analysis purposes we will not 
assume knowledge of the astrophysical mass distribution, 
and we will use a prior on the component masses that is 
uniform on the interval [1,2] Mq. As we shall see, sig¬ 
nificant biases will appear in the estimation of cq. These 
can be traced back to this flat prior. As demonstrated 
in the Appendix, if we had exact knowledge of the astro- 
physical mass distribution and could use that as a prior 
instead, the biases would go away. 


1. Zero spins; flat versus Gaussian distribution of 
component masses 

Let us start with the case of zero spins, and a uniform 
mass distribution. Fig. 8 shows the evolution of the me¬ 
dians and 95% confidence intervals in the measurement of 
Co as information from an increasing number of detected 
sources is combined, the injected EOS in turn being MSI, 
H4, and SQM3. We see that a clean separation between 
posterior densities occurs after ~ 50 sources have be¬ 
come available, and uncertainties of ~ 10% are reached 
as the number of detections goes towards 100. This can 
be compared with Fig. 1 of our earlier paper [25] . where 
the separation also happens around ^ 50 sources, but 
~ 10% errors are arrived at somewhat sooner than here. 
We recall that in that work, tidal effects were only taken 
to IPN order; on the other hand, waveforms were termi¬ 
nated at the LSO frequency rather than at the minimum 
of the LSO and contact frequencies. The earlier termina¬ 
tion of signal waveforms in the present paper leads to a 
smaller number of cycles, and somewhat less information 
about the EOS is available. 

In Fig. 9, we show results for zero spins, and this time a 
Gaussian distribution for the injected component masses. 
A good separation between MSI, H4, and SQM3 does not 
occur until ^ 150 sources have become available, and 
large systematic biases appear. As explained below, this 
is related to the continued use of a flat prior on the com¬ 
ponent masses, a distribution which now has a significant 
mismatch with the astrophysical one. The effect of the 
mass prior is further investigated in the Appendix. 



FIG. 8: Evolution of the medians and 95% confidence inter¬ 
vals in the measurement of co = A(mo), the tidal deformabil- 
ity at the reference mass mo = 1.35 Mq, for the cases where 
the injected EOS is MSI, H4, or SQM3. Both in the injec¬ 
tions and the templates, spins are set to zero, and the injected 
mass distribution is uniform on the interval [1, 2] Mq. 



FIG. 9: The same as in Fig. 8, but this time the signals have 
component masses drawn from a strongly peaked Gaussian 
distribution; on the other hand, the prior distribution for the 
masses used in the analysis of the data is still taken to be 
uniform on [1,2]M0. Note how large systematic errors ap¬ 
pear. The effect of the mass prior is further investigated in 
the Appendix. 


2. Gaussian mass distribution, non-zero spins 

We now focus on the case of a Gaussian distribution for 
the injected component masses, and also switch on spins, 
which are drawn from a Gaussian distribution with zero 
mean and a^ = 0.02. We also allow for spins in the tem¬ 
plate waveforms, with a prior distribution that is uniform 
on [—0.1,0.!], to reflect the ignorance of the true distri¬ 
bution of spins that we will have in reality. The results 
are shown in Fig. 10. As in the non-spinning case with 
the same injected mass distribution, there are system¬ 
atic biases. Having to estimate the spins as additional 
parameters also increases the statistical errors, because 
of the larger dimensionality of the parameter space to be 
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FIG. 10: The same as in Fig. 9, but now the signals not only 
have Gaussian distributed masses, but non-zero spins as well. 
Systematic errors remain, and statistical errors have increased 
due to the larger parameter space that needs to be probed. 

probed by the sampling algorithm. 

Finally, we mention that the higher-order coefficients 
Cl and C 2 are essentially unmeasurable in all the cases we 
considered (with or without a Gaussian mass distribution 
or spins); even with 100 sources, the posteriors are not 
significantly different from the priors. 


VI. DISCUSSION 

We have revisited the question of how well the equa¬ 
tion of state of neutron stars can be measured with obser¬ 
vations of binary neutron star inspirals using Advanced 
Virgo and Advanced LIGO. Our starting points were 
the Bayesian model selection and parameter estimation 
frameworks introduced in our earlier paper |25j . Given a 
set of hypotheses associated with a list of different EOSs 
one can calculate the odds ratios for all pairs in the set, 
which provides a ranking in which EOSs that are more 
similar to the underlying one will tend to come out near 
the top, whereas EOSs that differ from it signficantly will 
get deprecated. Another way to gain information about 
the EOS from multiple sources is to model the tidal de- 
formability A(m) as a series expansion in (m — mo)/MQ 
(with Too some reference mass), which is truncated at 
some suitable order. Since the coefficients in such an ex¬ 
pansion are source-independent, their posterior density 
distributions can be combined. For the EOS we con¬ 
sidered a “stiff” (MSI), “moderate” (H4), and “soft” 
(SQM3) equation of state, as well as the point particle 
model (PP). In it was found that for mg = 1.4 Mq, 
the deformability A(too) could be determined with ~ 10 % 
accuracy by combining information from 0 ( 20 ) sources. 
This was confirmed in recent work by Lackey and Wade 
[26j , who used a qualitatively similar waveform model as 
in |25j but implemented a more physical parametrization 
of the EOS in terms of piecewise polytropes. 

We have significantly extended our earlier study pS] . 


not only by expanding the number of simulated BNS 
sources, but also by incorporating as much of the rele¬ 
vant astrophysics as has been analytically modeled, such 
as tidal effects to the highest known order m , neutron 
star spins, the quadrupole-monopole interaction [laisD], 
the impact of possible early waveform termination due 
to the hnite radii of the neutron stars, and a strongly 
peaked Gaussian distribution of the component masses 

[SIHSll- 

In order to separate the impact of spins from the other 
effects, we first set spins to zero both in injections and 
templates (in which case the QM effect is also absent) 
while retaining the tidal effects as well as the potentially 
earlier termination of the waveform, and looked at hy¬ 
pothesis ranking for MSI injections. When choosing a 
wide, uniform distribution for the component masses, 
we saw that, as in [^, EOSs tend to be ordered cor¬ 
rectly according to stiffness and similarity to the true 
EOS. On the other hand, the log odds ratios between 
the incorrect and correct EOSs seemed to stretch to less 
negative values, presumably because of early waveform 
termination. Nevertheless (and again as in [IS]), hy¬ 
pothesis ranking worked well with catalogs of 0 ( 20 ) de¬ 
tected sources. The picture changed dramatically when 
the injected mass distribution was taken to be a strongly 
peaked Gaussian while keeping the mass prior to be uni¬ 
form and wide as before. In that case > 100 detections 
were needed to approach the discernibility of EOS seen 
in earlier work. Next we focused on a Gaussian distribu¬ 
tion for the masses, and switched on spins. At least for 
MSI injections, this turned out not to have a significant 
additional detrimental effect on our ability to distinguish 
between the EOSs. For H4, being in between MSI and 
SQM3 in terms of stiffness, we saw that the correct EOS 
got ranked above the others a reasonable fraction of the 
time, but the internal ordering became less clear. Finally, 
for SQM3, even with catalogs of 100 sources only MSI 
could be distinguished from the injected EOS reasonably 
well, but not H4 or PP. 

We also looked at parameter estimation for the coef¬ 
ficients in a series expansion of A(m) in the small quan¬ 
tity (to — too)/M 0 , truncated at some suitable order. 
Gontrary to our earlier work we used a quadratic rather 
than a linearized approximation; nevertheless we found 
that, here too, only the leading-order coefficient is mea¬ 
surable. When the signals have a strongly peaked Gaus¬ 
sian mass distribution rather than a flat one, again keep¬ 
ing the wide, flat mass prior, systematic errors are intro¬ 
duced. Switching on spins as additional parameters also 
increases the statistical errors. 

In the Appendix we investigated the effect on parame¬ 
ter estimation of the prior on the masses. We found that, 
if we can assume to have exact knowledge of the astro- 
physical distribution of the source masses so that it can 
be used as the prior distribution, the biases in the esti¬ 
mation of Co largely disappear. Recent estimates for this 
distribution [5Tti5i] are based on a rather small number 
of observed BNS systems and show dependence on the 
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methodology used; hence it seems that we cannot confi¬ 
dently claim to have detailed knowledge. One could con¬ 
sider supplementing the existing information with com¬ 
ponent mass measurements from the gravitational wave 
signals themselves, but as is well known, these will come 
with large uncertainties [33]; moreover, due to selection 
biases, the distribution of masses in electromagnetically 
observed neutron star binaries may differ from the mass 
distribution in BNS coalescences seen by gravitational 
wave detectors. A more extensive investigation of the ef¬ 
fect of the prior distribution of component masses is left 
for future work. 

There could be ways in which our conclusions are 
on the pessimistic side. For example, a more physical 
parametrization of the EOS as in allows one to fold 
in physical constraints such as causality, which is bound 
to improve parameter estimation. Moreover, it was re¬ 
cently found that the implementation of quantum squeez¬ 
ing in the interferometers may improve the measurability 
of tidal deformabilities by a few tens of percent |83|. Fi¬ 
nally, it is worth noting that with the “plausible” BNS 
detection rate of ~ 40 per year at design sensitivity [7], 
the desired number of sources could be collected over the 
course of a few years. 
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Appendix: Effect of the prior on component masses 


Unlike with our evidence calculations, in the case 
of posterior density functions it is relatively easy to 
“reweight” the sampling of parameter space so as to make 
p{0, {cj}\d, I) correspond to different priors on the pa¬ 
rameters [33]. The degradation in the estimation of cq 
( and for that matter, hypothesis ranking) happened when 
we changed the way the component masses in the injec¬ 
tions were distributed. Hence it is of interest to study 
the effect of the prior on the masses in particular. 


Let us pretend to have perfect knowledge of the astro- 
physical mass distribution - in our example a Gaussian 
with /im = 1.35 Mq and am — 0.05 Mq - and take the 
prior on mi, m 2 to be identical to it. In the case of zero 
spins, the result is shown in Fig. 11. We see that the sig¬ 
nificant biases we encountered in Fig. 9 have largely gone 
away. In Fig. 12 we also include spins as before; here too, 
the biases seen earlier are strongly mitigated, though the 
larger parameter space to be probed still causes larger 
statistical errors. 


This is not a typical case of a “prior-dominated” in¬ 
ference on a parameter, since the bias originates from a 
bad choice of priors for different parameters (mi, m 2 ) 
than the one that we are interested in (cq). Two im¬ 
portant details that make this bad choice manifest itself 
as a bias in the cq posteriors are the following. First, 
there is the fact that the parameters A^, through which 
Co is inferred, have an implicit dependence on the masses 
mA- The cq posterior is determined by the posterior on 
the m-A plane for each component NS, and if the masses 
are biased then so is the inferred A(m) curve. Second, 
since cq is treated as an independent parameter, the bias 
enters through the mass prior, in the process of marginal¬ 
izing over mi and m 2 , consistently for each source, and 
is therefore a persistent bias that will not average out as 
the number of sources increases. 


In conclusion, the biases we see in the estimation of 
Cq mostly result from the mismatch between the mass 
distribution for the sources and the prior distribution of 
component masses in the Bayesian analysis of the data. 
The relatively small remaining biases that occur when the 
injected mass distribution is the same as the prior can 
be attributed to the quadratic approximation for A(m) 
used in the template waveforms, and the fact that when 
most of the masses are in a narrow interval, less of the 
underlying tidal deformability function is being probed 
by the sources. 
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FIG. 11: The same as in Fig. 9, but this time using a Gaus¬ 
sian prior for the component masses that exactly matches the 
injected mass distribution. The significant biases that were 
seen before have largely disappeared. 
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Gaussian mass distribution 
r Gaussian mass prior 


Injected value 5QM3 
Injected value H4 
Injected value MSI 
95% Cl 5QM3 
95% Cl H4 
95% Cl MSI 



200 


FIG. 12: The same as in Fig. 11, but now with spins switched 
on. Again we use a Gaussian prior for the component masses 
that matches the injected distribution. Here too, the biases 
have been mitigated. 
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